function  [r12,r18,r1000]=SolveGeoSeries_shell(eta,theta)

%structural model
% eta=0.00512;
% theta= 0.05034;

beta=theta/eta;

 x = fsolve(@(x) x*(1-x^12)/(1-x)-beta,.1,optimset('Display','off'));
 
r12=1/x-1  ;
 
 x = fsolve(@(x) x*(1-x^18)/(1-x)-beta,.1,optimset('Display','off'));
 
r18=1/x-1  ;

 x = fsolve(@(x) x*(1-x^1000)/(1-x)-beta,.1,optimset('Display','off'));
 
r1000=1/x-1  ;

end